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ABSTRACT 

We have developed a simplified method of treating the radiative acceleration of dusty 
flows. This method retains the sharp impulse at the dust destruction radius that is a 
feature of frequency dependent radiative transfer, whilst placing minimal demands on 
computing resources. As such, it is suitable for inclusion in hydrodynamic codes. We 
have applied this method to the formation of massive stars in spherical geometry, and 
find that the fraction of a cloud which can accrete on to the central star is a strong 
function of the Jeans' Number and density profile of the cloud. Massive star formation 
is favoured by cold homogeneous conditions, as might result in regions where gas is 
swept up by some external triggering agent. We find (in contrast to previous authors) 
that massive star formation does not require a depleted dust mixture, although the use 
of dust at typical interstellar abundances does reduce the final stellar mass compared 
to cores formed from a depleted mixture. 

Key words: stars: formation - stars: early-type - hydrodynamics - methods: nu- 
merical - radiative transfer 



1 INTRODUCTION 

The question of massive star formation (> 10 M©) is of great 
astrophysical interest. Despite their rarity, massive stars 
dominate the light from young star clusters. As such, they 
are important tracers of star formation in distant galaxies. 
Within the Galaxy, massive stars also have profound effects 
on their local environment. During their brief lives, massive 
stars produce an intense UV radiation flux, ionis i ng an d 
dispersing nearby gas (see, e.g. Tenorio-Tagle et al, (198f:)), 
while powerful stellar winds are injected into the surround- 
ings. The death of a massive star is also a dramatic event 
- a supernova which disperses metals into the interstellar 
medium, while also injecting a large amount of mechanical 
energy. 

Given this wide range of consequences, it is evidently 
desirable to understand environmental effects that govern 
the formation of massive stars. Infall models for low mass 
stars have been studied for a number of years. Some of the 
earliest work in the area wa s perfo rmed (independently) by 
Larsorj ( |l969[ ) and [Penston] ( |196E| ). The possibility of mag- 
netic fields was considered by many authors, and ambipo- 
lar diffusion has been studied as a means of forming the 
quasi-static syst ems often assu med by star formation simu- 
l ation s (see, e.g. Nakano (1979) and Safier, McKee & Stabler 
(1997)). Infrared excesses of slightly older protostars (e.g. 



Mendoza (1968)) are generally interpre t ed as indicating the 
presence of accretion dis cs (see |Pringle| (|l98l[)). These ideas 
were drawn together by ^hu, Adams fc Lizanc (1987), and 
the infall model developed in this paper is now generally ac- 
cepted for low mass stars. However, massive stars require a 
more sophisticated treatment. 

The shorter evolutionary timescales of massive stars 
mean that they are likely to join the main sequence before 
the end of their main accretio n phase. This complicat e s sim - 



ulations of the core itself - see 



Tout, Livio fc Bonnell] ( ll99£| ) 
High 



for how this can be a problem even for low mass stars, 
accretion rates onto massive cores will also lead to large lu- 
minosities, which may significantly affect the accretion ffow 
through radiation pressure on dust grains. Moreover, since 
the main sequence lifetime of a massive star is comparable 
to the free fall time of the parent core, it is evident that the 
accretion history and luminosity evolution may be closely 
linked throughout the star's life. It should be noted that if 
stellar masses are self-limited by radiation pressure on dust 
grains, then the high mass end of the IMF is likely to be 
highly sensitive to the metallicity and dust grain properties 
of the parent star formi ng region. For more details, of these 



and other problems, see Stabler, Palla & Ho (2000) 
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The possibility of r adiation pressure ha l ting a ccretion 
was first consi dered by Larson fc Starrfield] (1971 ). Subse- 
quent work by [Wolfire fc Cassinellj ( [1986| , |1987| ) concluded 
that stars with masses > 100 M© could form, provided the 
dust abundance was depleted from normal galactic values. 
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Wolfire fc Cassinellil considered spherically symmetric ac- ties to compu te the radiative acceleration of the flow. Wolfire 



cretion onto a massive star and modeled the interaction be- 



Cassinelli 



tween the stellar/accretion luminosity and the infalling dust. 
Since interstellar dust typically melts when it reaches a tem- 
perature in the range 1500 — 2000 K, there is a significant 
amount of momentum transfer in a region close to the dust 
destruction radius, rd, where the radiation from the central 



core first impinges on the dusty inflow. Wolfire & Cassinelli 



used these models to propose several criteria necessary for 
the formation of massive stars, of which the two most im- 
portant are 

(i) that accretion should proceed at a sufficiently high 
rate that the ram pressure at rd exceeds the pressure asso- 
ciated with the radiation impulse at that point 

(ii) that the ratio of gravitational to radiative acceleration 
at the outer edge should exceed unity, in order that the 
initially static material will be set in motion in an inward 
direction 

These criteria have long been held to be inimical to the 
formation of massive stars and have prompted a range of 
alternative models, discussed below. The first criterion re- 
quires that massive stars form in dense cores (n ~ 10® cm~'^ 
for a 100 M0 core). This was regarded as uncomfortable 
because if such case was in a state of hydrostatic support 
prior to collapse (i.e. if the core contained roughly a Jeans' 
Mass) , the required temperatur e would be high (~ 1000 K 

which 



Wolfire & Cassinelli 



for a 100 M0 core according to 

is much higher than the temperature in such re gions - ac- 
cording to section 3.3 of Garay & Lizano (1999), tempera- 
tures are typically ~ 100 K, rising to a maximum of about 
250 K). The second criterion, on the other hand, translates 
(for a given stellar mass and luminosity) into an upper limit 
on the permitted value of the fiux mean opacity at the outer 
edge of the core. This turns out to be rather restrictive: for 
any reasonable value of the radiation temperature at the 
outer e dge, inflow is impossible for the s tandard dust mix- 
ture of Mathis, Rumpl & Nordsieck (1977). Accretion is only 



permitted if the dust mixture is substantially modified (re- 
duction in the total mass of grains by a factor of four from 
standard Galactic abundances - a factor of eight was used 
in their simulatio ns - and a smaller m aximum size of the 
graphite grains). Wolfire & Cassinelli thus concluded that 
'...special grain conditions must exist for massive star for- 
mation to occur.' 

These constraints have encouraged the discussion of al- 
ternative models of massive star formation. A straightfor- 
ward and realistic modification of the above is to allow for 
the angular momentum of the initial core and hence for the 
breaking of spherical symmetry due to the developme nt of a 



centri f ugally s upported disc. Yorke and col la borators ( Yorkc 
(tl979|, |1980|), [Bodcnhcimcr et all (ll990|), |Yorke, Boden 



heimer fc Laughlin ( 199^^ Yorke, Bodenheimer & Laughlin 



( 1995| ) and [Yorke fc Bodenheimer ( 199E )) have pioneered the 
hydrodynamic study of such a situation, including the effects 
of radiation pressure on dust. They have pointed out that 
the problem of massive star formation is considerably eased 
by the 'fiashlight' effect, whereby the bulk of the ultraviolet 
photons escape along the rotation axis of the disc, and thus 
do not intercept the accretion flow. For reasons of computa- 
tional economy, however, the majority of such calculations 
have been performed using grey (frequency averaged) opaci- 



have shown that this simplification may cause 
the efficacy of radiative feedback to be substantially under- 
estima ted (see below). Very recently, Yorke & Sonnhalter 
( [2002| ) have produced the first 2.5D hydrodynamic simula- 
tions in which the acceleration due to radiative feedback is 
computed using frequency dependent opacities. The com- 
plexities of such a code, and the large computational over- 
head it demands, has meant that it has not yet been possible 
to explore a large area of parameter space. Hence, it has not 
been possible to understand fully the relationship between 
initial conditions and final stellar mass. The continued in- 
crease in computational capacity should however bring this 
problem within gr asp in the next few years 



Alternatively, Bonnell, Bate & Zinneckei (1996) have 
suggested that massive stars form through a mixture of ac- 
cretion and coalescence in the cores of ultra-dense clusters. 
Coalescence is an attractive mechanism for massive star for- 
mation inasmuch as the effect of radiative feedback is neg- 
ligible. Nevertheless, in such simulations massive stars also 
acquire mass through Bondi-Hoyle type accretion and it is 
unclear how radiative feedback operates for this accretion 
geometry. 

In this paper, we revisit the issue of spherical accre- 
tion onto massive stars, through the use of hydrodynamical 
simulations in which the algorithm for radiative feedback 
is designed to match the outcome of frequency dependent 
radiative transfer calculations. We compare the re s ults o f 
such simulations with those of Wolfire & Cassinelli (1987), 



who constructed steady fiows on to stars of given mass and 
luminosity. We find significant differences in the results, in 
the sense that the conditions for massive star formation are 
apparently much le ss restrictive than those envisaged by 
Wolfire & Cassinelli: in particular, we find there is no need 



to modify the grain mixture, if collapse occurs from a highly 
Jeans unstable core. Needless to say, spherical accretion is 
not a realistic model for massive star formation, but it does 
represent the geometry in which radiative feedback is likely 
to be most effective (since both accretion and radiative feed- 
back occur over all solid angles). Since we have shown that 
the formation of massive stars by accretion is rather unprob- 
lematical even in spherical geometry, we may conclude from 
this that the obstacles to massive star formation are still less 
severe in the case of realistic (disc) geometries. 

The structure of this paper is as follows: Section ^ dis- 
cusses a simplified algorithm for radiative feedback, which 
preserves the features of a full frequency dependent treat- 
ment. In section ^we describe the hydrocode modified to in- 
clude the new algorithm. Results are presented in section ^ 
and a detailed discussion is in section |^. Conclusions and 
closing remarks are in section 



2 A NEW ALGORITHM FOR RADIATIVE 
FEEDBACK 

One of the most i mportant findings of the ra diative trans- 
fer calculations of Wolfire & Cassinelli (1987) was that the 



grey approximation is completely inadequate to describe the 
physics involved in the dusty infiow. The problem may be 



s een b y comparing figures 6 and 7 of Wolfire & Cassinelli 
(1987) - the calculations which used the Rosseland mean 
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opacity completely failed to produce the sharp deceleration 
at the inner edge of the dust shel l. This point is also made 
by Preibisch, Sonnhalter & Yorke (1995). 

The problem stems from the extreme manner in which 
the dust opacity varies with wavelength - typically k oc 
X~^'^. Use of the Rosseland mean opacity carries an im- 
plicit assumption of thermal equilibrium between the radi- 
ation and the fluid. At the inner edge of the dust shell, the 
dust (melting at ~ 2000 K) encounters the radiation field 
coming from the stellar surface (at a radiation temperature 
of ~ 2 X 10* K). The assumption of thermal equilibrium is 
evidently extremely poor at this point in the flow, and the 
steep variation of the dust opacity with wavelength ensures 
that this mismatch appea rs in the dynamics. Figure 8 of 
Wolfire & Cassinelli (1986) implies that the Rosseland op- 
tical depth of the envelope is no more than about 2, but it 
is easy to show that the optical depth at ultraviolet wave- 
lengths is on the order of hundreds. The use of Rosseland 
opacities therefore artificially softens the impact of the ra- 
diation field. 



2.1 The Simplified Model 

In order to remedy this situation, we here develop a simpli- 
fied algorithm which retains the sharp impulse at the inner 
edge of the dust shell (a feature of the full radiative transfer 
calculations), whilst at the same time effecting considerable 
computational economies. The prescription obtained is suit- 
able for incorporation into hydrodynamics codes. 



2.1.1 Radiative Transfer 



The e ssence of the method is to follow Wolfire & Cassinelli 
( [l986| ) in splitting the radiation field into two components: 
the direct stellar radiation field (L,) and a diffuse, ther- 
malised field (Lth). The total luminosity, comprised of the 
nuclear luminosity of the core and any accretion luminosity 
(which both change as the protostar evolves), is conserved, 
so that 

Ltot = L^, -\- Lth 



at all times and at all radii (when integrated over frequency) . 
The thermalised radiation field may be treated using grey 
opacities, since figures 6 and 7 of Wolfire & Cassinelli ( 1987 ) 
show that deep within the dust shell, this approximation is 
valid. 

The thermalised radiation field was calculated as fol- 
lows: At every radius, a pliotospheric temperature could be 
defined via 



(1) 



The radiation temperature thus obtained is that required 
to transport the radiative fiux, if the radiation became free 
streaming at that point. The optical depth to a potential 
photosphere at radius r may be calculated as 



r(r) = -KR(r,ad(r-)) / p(r') dr' 



(2) 



where, kr is the Rosseland opacity, and for simplicity, it is 
assumed that no significant reprocessing of radiation occurs 
outside the potential photosphere. The location of the ac- 
tual photosphere for the dust shell may then be found by 



the conventional definition: a photosphere lies at an optical 
depth of 2/3 from infinity. 

Outside the photosphere, the value of T^ad may be set 
equal to the radiation temperature of the photosphere, as 
per the assumption made by equation ^. Inward of the pho- 
tosphere, the diffusion approximation must be solved: 



dT,, 



dr 



I'th(r) 



(3) 



Since the grains are fiowing inwards, it is reasonable to as- 
sume that they will melt once the radiation temperature 
reaches Tsub ~ 2000 K, permitting the inner edge of the 
dust shell to be located (see below also). 

There remains the question of calculating Lth{r). This 
can be done by attenuating the stellar radiation field as it 
proceeds outwards from rd via 



Ll{r) ^ Ll{r^)e-^' 



where 



pAr 



(4) 



(5) 



is the wavelength-dependent optical depth of the inflow, and 
is the wa velength dependent r adiation pressure opacity. 



as defined by Wolfire & Cassinelli , There is the potential for 



this to have an unpleasant effect on equation ^, by affecting 
Trad - In this case, a slow iterative solution would be required. 
However, as outlined below, this problem was avoided in the 
simulations presented. 

The mechanical effect of the radiation field is provided 
using 



arad(r) = 



/feriA(r)dA 



(6) 



Anr'^c 

which may be calculated separately for the direct and ther- 
malised fields. Since the thermalised field is assumed to be a 
black body, the integral over wavelength may be calculated 
in a lookup table of Planck mean opacities. 

2.1.2 Dust Physics 

The dust physics were brutal l y pru ned, compared to the 
model of Wolfire & Cassinelli (1987). Instead of following 
the size evolution of the dust as it was destroyed, the dust 
is taken to be present or absent, based on its equilibrium 
temperature in the radiation field. This may be calculated 
using 



(47rrd) 



t{a,\)Ll dA = 



\a,\)Bx{T,^b)A\ 



(7) 

where Q^(a,A) is the absorption efficiency as a function of 
wavelength for grains of radius a. The desired sublimation 
temperature, Tsub, grain radius, and grain type may be se- 
lected appropriately. 

Melting the dust at a set temperature (and, implicitly, 
having it spontaneously reform at that temperature if out- 
fiow occurs) is a gross simplification to the true physical 
conditions. Howeve r , this a ssumption is often use d - see, e.g 



Bodenheimer et al, 
I 



nKruc 



( il990| ) ; [Preibisch etahl ( |l993| ) 



Bell et al 



has been calc ulated in studies of low mass stars 
Sedlmayr ( 199^ )) that, once con- 
ditions are favourable for the formation of dust, the grains 



(e.g. [Krueger, Ganger 
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grow r apidly. The caveat is the assertion of Sail & Sedlmayi 



The gas equation of motion was modified to be 



(1999) that nucleation is a rather difficult process. The ap- 



propriate destruction temperature is a matter of some de- 
bate. The temperature of 2000 K conventionally quoted is 
based on the sublimation of graphite grai ns. Ifowever, the 
dust destruction temperatures quoted by Lenzuni, Gail & 



Hennin 1995 ) are much lower - typically in the region of 



1200 K. Lenzuni et al, suggest that the discrepancy is caused 
by mista ken assumption s about dust destruction processes. 
However, Lenzuni et al. found that graphite grains tended to 
be destroyed by chemispluttering (collisions with hydrogen 
leading to the formation of short hydrocarbon chains) . Only 
silicates were destroyed by sublimation - and a t a na uch lower 
temperature. Duschl, Gail fc Tscharnuter (199C) present 
similar calculations for a protoplanetary disk. They also 
found that graphite grains were destroyed at relatively low 
temperatures by surface reactions, while the silicates subli- 
mated. However, their silicate grains survive to ~ 2100 K. In 
the dynamical collapses presented, Tsub = 1700 K was used. 

In the present work, all dust opacities were calculated 
from the efficien cy tables of Dr a ine. These tables are based 
on the work of 



(1993) and Laor 



Draine 



Draine & Malhotra 



Drainel (|l99j), and are updated versions 



of the values used by Wolfire & Cassinelli. The dust size dis- 



tribu tion wa s taken to be the standard power law of 
et al.| (H977|). 



MathiE 



2.1.3 Stellar Radiation 

In order to perform the radiative transfer calculation, it is 
necessary to specify the radiation incident on the inner edge 
of the dust shell. The spectrum of the stellar radiati o n fiel d 
was modified in the manner of [Wolfire fc CassineUj ( [l986| ), 
with radiation more energetic than the Lyman Limit repro- 
cessed to lie below it. That is, the stellar spectrum, BX is 
chosen such that 



Bl 





eBx 



for < A < 912 A 
for A > 912 A 



and 



/ BxA\= / 
Jo Jo 



Bl dA 



(8) 



(9) 



where B\ is a true black body spectrum, and e is chosen to 
ensure equation |^ remains true. Then, 



LI = 47v^BlRl 



2.2 Testing the Simplifications 

The approximations described above w ere tested against the 
full radiative transfer calculations of 

The reference model was the 100 Mr calculation of 



([1987D 



Wolfire & Cassinelli 



Wolfire 



CassineUi (1987). This model had a steady ac- 
cretion rate of 5 x 10"'^ Mq yr"^ and a total core luminos- 
ity of 2.43 X 10® Lq. Where possible, boundary conditions 



were taken to be those specified by Wolfire & Cassinell 
so equation |^ was not used. Use of these boundary condl 
tions also eliminated the troublesome relation between equa- 
tions ^ and ^, since the dust would always be destroyed in 
the innermost grid cell. 



1 du^ 
2"d7 



GM{r) IdP Jkl'Lxir)dX 



p Ar 



+ 



(10) 



where M(r) is the mass enclosed at radius r, and P is the gas 
pressure. The A subscripts indicate quantities that ar e wave- 
length dependent. This is based on equation 31 of Wolfire 



CassineUi (1987). However, the left hand side has been 



changed into a ram pressure gradient, while the dust to gas 
mass ratio (which is ~ lO"'') has been dropped from the 
right hand side. 

The steady state velocity field obtained is plotted in 
Figure ^. It may be seen that, although the general be- 
haviour has been duplicated, there is a discrepancy of a 
few kilometres per second in the inner portions of the flow. 
This difference can probably be ascribed to the extremely 
approximate nature of the radiative transfer method used. 
Since the dust shell has an optical depth of the order unity at 
temperatures appropriate to the diffuse field, it is inevitable 
that errors will creep in - this regime is notoriously tricky to 
approximate. While it is certain that some method could 
be concocted to improve this specific case, its generality 
would be suspect. It should be noted that the magnitude 
of the jump in the infiow veloc ity at the inner edge of th e 
dust shell is identical to that of Wolfire & CassineUi (1987). 
Comparing the momentum fluxes of the curves, the initial 
encounter with the stellar radiation gives a change equiv- 
alent to 0.96Z//c, while the thermalised field adds another 
Q.7\L/c. These numbers are as one might expect - all the 
luminosity^ is absorbed and re-emitted at the inner edge of 
the dust shell. However, the thermal processing is signifl- 
cant, contributing almost as much momentum fiux as the 
initial strike. Again, this is as one would expect, given that 
the entire shell has tr ~ 1. 

The algorithm presented is a substantial improvement 
over calculations using only the Rosseland mean (which can- 
not even reproduce the qualitative behaviour of the inflow). 
Although there is a discrepancy of a few kms"^ at the inner 
edge of the dust shell, we note that the regime under con- 
sideration was the most problematic to approximate, with 
Rosseland optical depths only slightly larger than unity. We 
would expect the approximations we have made to become 
better as the dust shell becomes optically thicker. 



3 THE DYNAMICS CODE 



Stone & Norman 



( [1992D was cho- 



The ZEUS-2D code of 

sen to simulate the collapses. Since the flows were generally 
found to be highly supersonic, and fairly optically thin in 
the Rosseland mean, the gas equation of state was taken to 
be isothermal, with the temperature set at the start of the 
simulation. 

Resolution considerations lead to a slight modification 
of the algorithm detailed in section ^ The loop between 
equations ^ and ^ was broken by noting that the fiow is 
likely to be so optically thick that the direct stellar field will 
be attenuated in a single grid cell. In this case, Lt\i{r) — Ltot 

^ The missing four percent is probably due to the changing ve- 
locity affecting the density structure (since the accretion rate is 
constant), which will then change the gravitational forces slightly 
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Distance from centre /c 



Figure 1. Inflow velocities onto 100 Mq core. The thick soli d curve is that obtained fro m the new algorithm, while the thin solid curve is 



the freefall solution. The dashed curves are points read from 



Wblfire & CassincUi 



(1987), with the short dashes following the full solution 



of their figure 6, while the long dashes follow the Rosseland mean calculation of their figure 7. The dotted line is the curve one obtains 
if the effects of Lth are omitted in the new algorithm 



outside this grid cell. The mechanical impulse exerted on the 
inflow at the dust destruction radius may be represented by 
the injection of Ltot/c of momentum into the dust destruc- 
tion grid cell. If a photosphere is not found according to 
equation 2, then the value of may be calculated using 
equation i and the stellar field swept outwards from this 
point using equation |^. No thermalised field need be calcu- 
lated in this case, since its effect will be minimal. 

To avoid the central grid cell reaching unreasonable den- 
sities (and thereby consuming impractical amounts of com- 
puter time), ZEUS was modified so that a linearly decreas- 
ing amount of mass was accreted from the innermost grid 
cells, and placed into a point mass at the centre of the grid. 
This simulated the formation of the protostar - which is not 
resolved in our work. The number of cells which were per- 
mitted to accrete, and the maximum fraction accreted were 
found to have little effect on the results. The point mass 
was used to determine the intrinsic luminosity and the ra- 
dius (for the accretion luminosity) of the forming protostar. 
These quantities were determi ned using the Zero Age Main 
Sequence (ZAMS) formulae of [Tout et all ( |l996[ ).^ To allow 
for the fact that the core would not imr nediate l y sett le to 
its main sequence radius, the method of Vorke (197£) was 
adopted: once the core reached 0.1 Mq, its radius was set 
equal to that of the first grid cell. It was then allowed to 
contract on its instantaneous Kelvin-Helmholtz timescale, 



until it shrank to its main sequence radius. It is of course 
undesirable for the physical behaviour to be determined by 
computational considerations (the positioning of the grid 
cell). However, tests at different resolutions found no signif- 
icant variation in behaviour due to this. 

Most runs were performed with a grid of 600 cells, fifty 
of which were concentrated in the inner 100 an with linear 
spacing. The remainder were logarithmically spaced. Higher 
resolution runs were also performed, typically with 800 grid 
cells following the same general distribution pattern. 



4 RESULTS 

Two sets of simulations were performed. The first started 
from uniform initial conditions, the second started from a 
power law density profile. 



4.1 Collapse from Uniform Initial Conditions 

In Tables g, I and I we detail a number of simulations which 
start from uniform initial density. The initial density was 
10"^" gcm""^, with the different masses obtained by varying 
the simulation volume. The temperature of the gas, Tgas of 
each run was fixed at the indicated temperature; the corre- 
sponding value of a is also listed, where a was defined as 



^ Note that this assumption will overestimate the luminosity of 
the core, and increase the effects of radiative feedback 



Eth 



(11) 
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Tgas/K 


a 


M./Mq 


Comments 


Tgas/K 


a 


M,/Mq 


Comments 


10 


0.03 


196 


No thermalised field 


10 


0.05 


64.0 




200 


0.54 


88 


No thermalised field 


20 


0.11 


56.0 




10 


0.03 


185.7 




50 


0.27 


40.5 




10 


0.03 


188.3 


High resolution 


100 


0.54 


30.2 




10 


0.03 


183.8 


High resolution 


150 


0.81 


26.5 




20 


0.06 


172.6 




200 


1.08 


26.8 




50 


0.14 


126.5 




300 


1.62 




Did not collapse 


100 


0.27 


84.4 




10 


0.05 


37.6 


Normal dust levels 


150 


0.41 


67.1 




20 


0.11 


23.7 


Normal dust levels 


200 


0.54 


48.7 




50 


0.27 


12.2 


Normal dust levels 


200 


0.54 


47.7 


High resolution 


100 


0.54 


11.5 


Normal dust levels 


300 


0.81 


39.9 




150 


0.81 


14.2 


Normal dust levels 


400 


1.08 


45.1 




200 


1.08 


11.8 


Normal dust levels 


10 


0.03 


137.1 


Normal dust levels 


300 


1.62 




Normal dust levels. 


10 


0.03 


138.4 


Normal dust levels, 








did not collapse 


20 


0.06 


106.4 


high resolution 
Normal dust levels 


Table 2. Results from 88 Mq collapses of uniform gas 


50 


0.14 


54.5 


Normal dust levels 










100 
150 
200 
300 
400 


0.27 
0.41 
0.54 
0.81 
1.08 


15.8 
13.6 
12.7 
12.3 
12.4 


Normal dust levels 
Normal dust levels 
Normal dust levels 
Normal dust levels 
Normal dust levels 


Tgas/K 


a 


Af*/M0 


Comments 


10 
200 

10 
200 


0.003 
0.067 
0.003 
0.067 


4809 
3976 
4368 
2373 


No thermalised field 
No thermalised field 


400 


1.08 


11.7 


Normal dust levels, 










high resolution 


Table 3. Results from 5630 M© collapses of uniform gas 



Table 1. Results from 241 Mq collapses of uniform gas 



where Sgrav is the system gravitational energy and 

37^A/^gas 

where TZ is the molar gas constant.^ Note that by setting 
cv = §7?. the freeze-out of the rotational modes of molec- 
ular hydrogen is ignored. However, this is appropriate to 
the isothermal nature of the simulations. Unless stated oth- 
erwise, the dust mix t ure is the depleted one employed by 
|l987). The 'normal' dust mixture is 



Wolfire fc Cassinelli 

that of [ Praine fc Lee| ( |1984| ), which is used by [Wolfire & 
Cassine^lil ( |1986| ). Some of the simulations had no ther- 
malised radiation field present (that is, Lth = everywhere), 
as indicated. It should be noted that these runs lead to 
higher final stellar masses than their counterparts with the 
full feedback mechanism, pointing to the importance of the 
reprocessed radiation. Variation of the grid resolution was 
not found to have a significant impact on the final stellar 
masses. 

Figures ^ ^, ^ and |^ depict snapshots of the flow pro- 
file and illustrate the chief phases of the collapse of homo- 
geneous clouds subject to radiative feedback. In all cases, a 
core forms after a free-fall time, where 



behaviour discussed by Larson (196!:) and |^enston (196£) 
- a flat central density proflle with r^^ wings (Figure 
The colder collapses produce larger constant density regions 
£18 they approach homologous collapse, with steeper wings 
(typically closer to , see Figure 

Once the core has formed, the inner region of the flow, 
which is dominated by the core's gravity, approaches the 
r^^''" profile expected for steady flow onto a point mass (Fig- 
ure ^). At this stage, the dominant contribution to the core 
luminosity is initially from accretion, but the relative impor- 
tance of the intrinsic core luminosity increases as the core 
mass grows (Figure ^). In the flnal phase of the evolution, 
the accretion flow starts to be modulated by the effect of 
radiative feedback on the dust - oscillations in the accretion 
rate may be seen in Figure ^ Radiation pressure temporarily 
reduces the flow rate from the dust destruction radius, which 
produces a drop in the accretion rate onto the core after one 
freefall time from rd. The corresponding drop in luminosity 



32Gpo 



= 2.1 X 10' 



10" 



gem 



which is the standard result for the pressure free collapse 
of a homogeneous cloud. There is obviously no core lumi- 
nosity (or feedback applied) up to the point of core forma- 
tion. The density profile during the initial collapse depends 
on the temperature of the flow: warmer flows exhibit the 

^ Note that with this definition, a is related to the Jeans' Num- 
_ 3 

ber, Nj hj Nj = a 2 




Figure 2. Initial collapse of 200 K 241 Mq cloud. The curves are 
plotted for t = and after 0.17, 0.34, 0.51, 0.69, 0.84, and 1.0 
freefall times. Dotted line is curve 
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Figure 3. Initial collapse of 10 K 241 Mq cloud. The curves are 
plotted for t = and after 0.76, 0.84, 0.94 and 1.0 freefall times. 
Dotted line is r~'^ curve 



Figure 6. Accretion rates for 10 K and 200 K 241 Mq clouds. 
Solid line is for 10 K gas, dotted line follows 200 K cloud 




Figure 4. Middle phase of collapse of 10 K 241 Mq cloud. The 
core grew from 101 Mq to 156 Mq during the period of 381 yr 
plotted. Dotted line is r~^-^ curve 



causes a reduction in feedback at rd, and the accretion rate 
climbs again. The effect of this time delayed feedback is im- 
printed as oscillations in the density profile (Figure |^. In 
the last profile, the large drop in density at the inner edge 
of the dust shell is plainly visible, as the flow is effectively 
stalled by the momentum imparted by the photon field. 
At a quantitative level, the behaviour of the fiow can 




Figure 5. Comparison of intrinsic (nuclear) luminosity to to- 
tal core luminosity for 10 K (solid line) and 200 K (dotted line) 
241 Mq clouds 



Figure 7. Terminal phase of collapse of 10 K 241 Mq cloud. This 
plot covers a period of 1200 yr at the end of the simulation 



be understood in terms of the analytic expressions presented 
in Appendix ^ In particular, accretion is found to stall at 
the point where the momentum fiux in the accretion flow 
at rd matches L/c from the core. The evolution of the dust 
destruction radius (computed in the code according to the 
algorithm of section ^ may be compared with the approx- 
imation L = Ar^ made in the Appendix, and was found 
to be in satisfactory agreement over most of the evolution. 
Values of A were found to be similar to those of Wolfire &i 



Cassinelli 



198' 



As Figure H demonstrates, the colder gas clouds un- 
derwent a spike of accretion as shells from a range of radii 
arrive almost simultaneously at the centre (as expected for a 
homogeneous collapse). Roughly half the mass of the cloud 
goes into the core during this first spike. The associated ac- 
cretion rate is high, and briefiy super-Eddington - see Fig- 
ure For a star of mass M Mq, the Eddington Luminosity 
is 3 X lO^Af Lq. After the initial accretion spike, accretion 
continued but soon large oscillations develop in the accretion 
rate, and accretion subsequently stalls. 

The 'warm' run behaves somewhat differently. The ac- 

^ Note that our simulation takes no account of electron scat- 
tering opacity, since the ionised regions of the flow arc not re- 
solved. However, since the integrated energy input during the 
super-Eddington phase is rather less than the binding energy of 
the core, the core would not be liable to disruption 
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Figure 8. Comparison of total luminosity to Eddington luminos- 
ity for 10 K (solid line) and 200 K (dotted line) 241 Mq clouds 



Figure 9. Accretion histories for cloud started from initial 
density profile and no feedback. Solid curve is for 10 K gas, the 
dotted line 100 K, and the dashed line 200 K gas 





a 


M,/Mq 


Comments 


10 


0.02 


25.3 




20 


0.05 


28.0 




50 


0.11 


28.2 




100 


0.23 


36.3 




150 


0.35 


51.4 




200 


0.46 


64.8 




300 


0.68 


53.2 




400 


0.92 


45.7 




10 


0.02 


14.3 


Normal dust levels 


20 


0.05 


15.8 


Normal dust levels 


50 


0.11 


15.4 


Normal dust levels 


100 


0.23 


11.5 


Normal dust levels 


150 


0.35 


12.0 


Normal dust levels 


200 


0.46 


11.3 


Normal dust levels 


300 


0.68 


11.5 


Normal dust levels 


400 


0.92 


11.7 


Normal dust levels 



Table 4. Results from 262 Mq collapses of p oc r ^ gas 

cretion rate rises more slowly, and remains steady for about 
10% of the freefall time of the original cloud. The rising 
intrinsic luminosity of the core causes feedback to play an 
increasing role. Oscillations develop in the accretion rate, 
leading to the eventual stalling of the accretion flow. 

4.2 Collapse from Power Law Initial Conditions 

We now consider collapses in which the initial density profile 
scales as . Final stellar masses as a function of a (still 
relevant for this density distribution) are given in Table ^ 
The densities were normalised so that models of given mass 
and temperature had a values similar to those starting from 
uniform conditions. 

The effect of feedback may be most readily understood 
by first considering the growth of the central core in such 
models when feedback is not included. A set of sample curves 
is presented in Figure ^. The core growth of the coldest run 
roughly follows the M oc law derived by considering the 
freefall collapse of thin shells. The core is therefore built 
over a much longer interval than the homogeneous case. The 
warmer runs behave in a manner similar to their homoge- 
neous counterparts, since early core formation is suppressed 
by the action of pressure gradients. Since M{r) oc , the 
gravitational acceleration, g, is constant throughout the ini- 
tial cloud. However, the acceleration due to the pressure 




Figure 10. Collapse of 200 K cloud started from initial den- 
sity profile. The solid curves are plotted for 0.36, 0.72, and 0.86 
freefall times from the edge of the simulation. The dashed line is 
the initial density profile. Dotted line is curve 



gradient varies according to 

p dr p r r 

so pressure support improves towards the centre of the cloud. 
Consequently, the gas evolved towards the Larson-Penston 
solution initially, as maybe seen from Figure The subse- 
quent behaviour of the warmer clouds closely mirrored those 
collapsing from uniform initial conditions. 

When feedback was included, it was found that the 
warmest runs led to final stellar masses almost indistinguish- 
able from the corresponding homogeneous simulations (c.f. 
Tables |l| and ^). This is as one would expect, given the ap- 
pearance of the Larson-Penston solution prior to core for- 
mation. However, for still lower temperatures (Tgas < 200 K, 
or a < 0.5 for the conditions considered here), the trend is 
reversed, and final stellar mass decreases with lower temper- 
atures. This behaviour is explained by Figure ^. The lack of 
pressure support caused the cores to form over a longer pe- 
riod of time, since the formation of the core is not delayed by 
pressure gradients. This permits feedback to operate when a 
substantial fraction of the initial cloud mass still lies beyond 

The interaction of the collapse with the radiative feed- 
back followed a similar pattern to that observed in the col- 
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lapses from uniform clouds. The collapses with feedback 
closely followed the behaviour of those without, until the 
central luminosity became significant. Then, the flow be- 
gan to oscillate on a short timescale, and accretion halted 
rapidly. 



5 DISCUSSION 

Three phases of collapse were observed, where a different 
physical effect was dominant: 



(i) Self-gravity of the cloud 

(ii) Gravity of the core 

(iii) Luminosity of the core 



Stages |(i)j and 'ii) have been studied for a number of years. 
Radiative feedback is important in the final phase. 

A cursory glance at Tables 0, ^ and ^ suggests that, 
given the right initial conditions, massive stars may form 
despite the effects of radiative feedback. For example. Ta- 
ble M documents the formation of stars containing thousands 
of solar masses of material; we cannot of course rule out the 
possibility that such objects would fragment into a cluster 
of lower mass stars, but this is not a radiative effect. Such 
stars are not observed, but we find no reason why radiative 
feedback should prevent their formation. Higher initial cloud 
masses lead to higher final stellar masses. However, we do 
find that values of a (defined in equation ^) close to unity 
cause the fraction of the cloud mass accreted to be very low, 
whereas small a values gave rise to a high accretion frac- 
tion. This variation with a ar ises because feedback can only 
become effective in stage (iii) of the collapse. For the colder 



homogeneous clouds, most of the mass had already arrived 
within rd by the time this stage was reached. By contrast, 
warmer collapses possessed pressure gradients which mod- 
ified the density profile, and caused the core to form over 
a more extended period. Indeed, the steady accretion rate 
of the warm cloud depicted in Figure ^ is rather simila r 
to the conditions considered by [Wolfire fc Cassinellil ( [l987| ), 
although our differing parameterisations of the protostellar 
core prevent quantitative comparison. By contrast, the cold 
collapse is qualitatively different and points to the possibil- 
ity that massive star formation may not be as difficult as 
previously envisaged 



' Woltire & (Jassmelli (1987) concluded that depleted 
dust abundances were necessary to form massive stars, yet 
the current work has demonstrat e d the formation of stars 
> 100 M0 for the Draine & Lee ( 1984|) dust abundanc es. 

to 



Wolfire & Cassinelli 



However, the condition derived by 
support this requirement was analogous to the Eddington 
Limit, and hence was based on accelerations. Balancing ac- 
celerations is only useful when the processes involved are 
close to equilibrium (e.g. radial motion through an accre- 
tion disc). The present work studied highly dynamical col- 
lapses, and the velocities obtained in the early stages of the 
collapse were easily sufficient to overcome the net outward 
acceleration as the luminosity rose. 

A plot of initial a against the fraction of mass accreted 
onto the star (a measure of the efficiency) for an initially 
uniform cloud is presented in Figure |ll|. In general, lower 
Q values lead to more massive stars. For q ~ 1 and higher 
dust opacities, the final stellar mass obtained seems to tend 



to a constant value (c.f. Tables ^ and ^), possibly indicating 
that a true radiation pressure limit is being found. For the 
depleted dust mixture, the stars obtained can still be quite 
massive - > 40 M© . 

The results from collapses of p oc gas, as presented 
in Table ^ show that radiative feedback can be very im- 
portant under more condensed initial conditions. A plot of 
a against the star formation efficiency is presented in Fig- 
ure 12. The rise in efficiency with increasing temperature 
is p ain for the depleted dust, as is the steady value ob- 
tained for normal dust abundances. In the case of normal 
dust abundances, a limit of ~ 12 Mq is implied, with the 
variations around this value not being very significant. For 
the depleted dust mixture, the rather surprising result that 
colder gas does not lead to more massive stars is related 
to the behaviour of the collapse before radiative feedback 
becomes significant. If thermal effects are significant, they 
can delay core formation, and hence the onset of radiative 
feedback. During this delay, gas settled inside the dust de- 
struction radius (see Figure M and associated discussion) . 



6 CONCLUSION 

We have developed a new simplified algorithm for calculat- 
ing radiative feedback in hydrodynamic simulations. This 
method combines the computational economy of a frequency 
averaged treatment of the thermalised field (via a modified 
diffusion approximation) with the benefits of a frequency de- 
pendent treatment of the attenuated stellar field. This lat- 
ter feature ensures that the flow receives the sharp impulse 
at the dust destruction radius that is found in full radia- 
tive transfer calculations; we have demonstrated (Section ^ 
that this algorithm is significantly superior to methods us- 
ing grey opacities, which have become the standard method 
in hydrodynamical simulations to date. 

We have applied this algorithm to the formation of mas- 
sive stars by spherically symmetric accretion. Although this 
obviously does not represent a realistic geometry for OB star 
formation, our motivation has been to compare the hydro- 
dynamical trea^msnt_of_tlu^^ the steady state 



models of Wolfire & Cassinelli (1987), in order to evaluate 



how the time-dependence of the hydrodynamic calculation 
affects the results. Since we find (see below) that massive 
star formation is rather unproblematical (given appropriate 
initial conditions), we would expect it to be yet easier with 
a realistic (disc) accretion geometry. Our chief conclusions 
are as follows: 

(i) Massive star formation is favoured by models in which 
the luminous core forms quickly. In this case, much of the 
cloud is already within the dust destruction radius (and 
hence immune from feedback) when the core luminosity be- 
comes significant. The most favourable conditions involve 
cold homogeneous collapse, since in this case the core as- 
sembly time is shortest. Warmer collapses (with thermal to 
gravitational energy ratios near unity) produce lower final 
stellar masses, but can still readily produce OB stars given 
sufficiently massive progenitor clouds. For these warm col- 
lapses, the final stellar masses are rather insensitive to the 
initial cloud density profile 



(ii) We find that, in contrast to Wolfire fc Cassinelli 
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Figure 11. Initial a against final fraction of mass accreted for collapses from a uniform cloud. Note that Nj R! a 2 for the uniform 

collapses 
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Figure 12. Initial a against final fraction of mass accreted for collapses from an initial density profile of p ex r ^. In all cases, the total 
cloud mass was 262 Mq 
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(1987), a depleted dust mixture is not necessary for the for- 
mation of OB stars. The reason for this difference is that 



which leads to 



Wolfire & Cassinelli required that the net acceleration at 



the cloud's outward edge was inward, in order to attract the 
initially stationary material. In the hydrodynamic collapses, 
the fluid is already inflowing at close to its freefall velocity by 
the time the core luminosity becomes significant. However, 
increasing the optical depth implies a stronger deceleration 
of the flow from the difi'use field, which does reduce the final 
stellar masses. 

Although, as stated above, we do not consider our sim- 
ulations to be a realistic portrait of massive star formation, 
we end with a few remarks about the possible relevance of 
these findings to OB star formation. Firstly, we might expect 
that the qualitative conclusion that massive star formation is 
favoured by rapid core formation to hold in other geometries. 
Such rapid formation is promoted by rather homogeneous, 
highly Jeans unstable conditions - as might be expe cted to 
arise in regi ons of triggered star formation (see, e.g. [Braun 
" (|l997| ) 



Olsen, Kim fc Buss (2001)). Secondly 



et al 

note that in our models in which feedback effects were par- 
ticularly severe (either due to warm initial conditions or the 
use of a standard grain opacity), the resulting stellar mass 
often ended up near to ~ 12 Mq. The lack of any obvious 
signature in the IMF around this value therefore argues that 
the dominant population of progenitor clouds does not fall 
into this category. 



APPENDIX A: SCALING RELATIONS 

For accretion to occur, despite radiation pressure on the em- 
bedded dust grains, the ram pressure at the dust destruction 
radius, rd must exceed the momentum of the radiation field. 
That is 



L 

Mui > — 
c 



(Al) 



where tt d is the inflow veloc i ty at rd. This is the condition 
used by [Wolfire fc CassineUi ( 1987 ) at the inner edge of their 
dust shell. 

The value of Ud may be calculated on the assumption 
that the stellar mass dominates, and that the material is in 
free fall at the destruction radius. To make further progress, 
assumptions must be made about the source of the luminos- 
ity. 



as the condition for accretion. That is, the dust destruction 
radius must be less than some increasing function of the 
stellar mass. This would mean that accretion becomes eas- 
ier over time, and that radiative feedback would place no 
limit on stellar masses. Evaluating the constants gives the 
following prediction for steady state accretion to occur: 



rd < 4400 



(A2) 



Of course, the dust destruction radius, r^ varies with 
core luminosity. As a first approximation, the dust destruc- 
tion radius will vary as 



Arl 



(A3) 



for some constant A. Table 3 of Wolfire & Cassinelli (11987 

. _ O 1 ' J . . ' 



suggests that A « 3 x 10* ergs"^ cm . If it is assumed that 



equation 



A3 is an expression of the Stefan-Boltzmann Law, 

far below 



then the implied temperature is less than 1000 K ■ 
the nominal dust melting temperature, and another exam- 
ple of the importance of the wavelength dependence of the 
dust opacity. However, with this dependence folded into A, 
equation A3 may be safely used to estimate the value of rd. 

Putting this together with the free fall assumption 
yields 

Li 



MV2GM, > 



Ad^ 



(A4) 



Substituting an accretion luminosity then yields 

M < 'iG''-^M:^^RlAc' 
Finally, putting in the mass-radius relation gives 
AARq'^c'^ 



M < 



G3Mn3 



6 Mq yr-i 



(A5) 



as the condition for accretion. Rather surprisingly, this is 
independent of the stellar mass. However, no star of reason- 
able mass could accrete at this rate and remain on the main 
sequence. We therefore conclude that accretion luminosity 
on its own is unlikely to shut down an accretion flow. 



Al Luminosity Dominated by Accretion 

If the core's luminosity is dominated by accretion, one ob- 
tains 

MV2GM. > — r 

cR, 

To compute the accretion luminosity, a mass-radius re- 
lation is required. A simple mass radius relation for massive 
stars is 



A2 Luminosity Dominated by Nuclear Burning 

If the luminosity is dominated by nuclear reactions in the 
core, a different relation between L and M is obtained. For 
main sequence stars, L ~ , where ip is in the range of two 
to three for most stars (although the most massive stars have 
V' ~ 1). Substituting the power law relation into equation A4 
suggests that the condition on the accretion rate will become 



^0 



4G2ylc4Mo5V' 



M. 



5i>-2 



(A6) 



Within the likely range of values for t/i, this minimum ac- 
cretion rate is a rather rapidly increasing function of M». 
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Hence, accretion is becoming harder, and is likely to shut 
down. 
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